// vim:et:ts=2:sw=2:ci:cino=f0,g0,t0,+0:

#define START 0

#include "TMV_Test.h"
#include "TMV_Test2.h"
#include "TMV.h"
#include "TMV_Sym.h"
#include "TMV_TestSymArith.h"

#define NOLDIVEQ
#define NORDIVEQ
#include "TMV_TestMatrixDivArith.h"

template <class T> void TestSymDiv_A(tmv::DivType dt, PosDefCode pdc)
{
  const int N = 10;

  std::vector<tmv::SymMatrixView<T> > s;
  std::vector<tmv::SymMatrixView<std::complex<T> > > cs;
  std::vector<tmv::BaseMatrix<T>*> B;
  std::vector<tmv::BaseMatrix<std::complex<T> >*> CB;
  MakeSymList(s,cs,B,CB,pdc);

  tmv::Matrix<T> a1(N,N);
  for (int i=0; i<N; ++i) for (int j=0; j<N; ++j) a1(i,j) = T(1-3*i+j);
  a1.diag().AddToAll(T(10)*N);
  a1 /= T(10);
  tmv::Matrix<std::complex<T> > ca1 = a1 * std::complex<T>(3,-4);

  tmv::SymMatrix<T> s1(a1);
  tmv::SymMatrix<std::complex<T> > cs1(ca1);
  tmv::HermMatrix<T> h1(a1);
  tmv::HermMatrix<std::complex<T> > ch1(ca1);

  tmv::SymMatrix<T> s1x = s1;
  tmv::SymMatrix<std::complex<T> > cs1x = s1;
  tmv::HermMatrix<T> h1x = h1;
  tmv::HermMatrix<std::complex<T> > ch1x = ch1;

  for(size_t i=START;i<s.size();i++) {
    if (showstartdone)
      std::cout<<"Start loop: i = "<<i<<", si = "<<tmv::TypeText(s[i])<<"  "<<s[i]<<std::endl;

    tmv::SymMatrixView<T> si = s[i];
    tmv::SymMatrixView<std::complex<T> > csi = cs[i];
    if (dt == tmv::CH && csi.issym()) continue;
    si.SaveDiv();
    csi.SaveDiv();


    TestMatrixDivArith2<T>(dt,s1x,cs1x,si,s1.View(),csi,cs1.View(),"Sym/Sym");
    TestMatrixDivArith1<T>(dt,h1x,ch1x,si,h1.View(),csi,ch1.View(),"Herm/Sym");
  }
  for(size_t i=0;i<B.size();++i) delete B[i];
  for(size_t i=0;i<CB.size();++i) delete CB[i];
}

#ifdef INST_DOUBLE
template void TestSymDiv_A<double>(tmv::DivType dt, PosDefCode pc);
#endif
#ifdef INST_FLOAT
template void TestSymDiv_A<float>(tmv::DivType dt, PosDefCode pc);
#endif
#ifdef INST_LONGDOUBLE
template void TestSymDiv_A<long double>(tmv::DivType dt, PosDefCode pc);
#endif
